program PS

        ! "use" is the inheritance mechanism for F90:
        use mConstants
        use mObject
        use mCosmology
        use mLPS
        use mPSFittingFunctions
        use mMisc

        implicit none

        real(D), parameter :: K_LL = 0.01_D, &
                              K_UL = 1000.0_D, &
                              SCALE_FACTOR = 0.5_D

        integer, parameter :: NUMPTS = 100

        real(D) :: k(NUMPTS), lightSmithPS, darkSmithPS
        integer :: i

        ! "type" is equivalent to a C struct:
        type(tUniverse) :: darkUniverse
        type(tLPS) :: darkLPS
        type(tSmith) :: lightSmith, darkSmith

        ! Think of "new" as a constructor. This makes a LightSmith "object" using all the default cosmology parameters in base/object.f90:
        call new(lightSmith)

        ! If you want to change some parameters you can do it like this:
        call new(darkUniverse, w_0 = -0.8_D, w_a = -0.6) ! A new universe with dark energy.
        call new(darkLPS, universe = darkUniverse, normalization = 'CB') ! Linear power spectrum with dark energy and CMB normalization.
        call new(darkSmith, darkLPS)

        write( 6, '(a)' ) '# k(h/Mpc) P_Smith P_Smith_DE'

        call geometricPartition( k, K_LL, K_UL )

        do i = 1, NUMPTS

                call smithEtAl( lightSmith, k(i), SCALE_FACTOR, lightSmithPS )
                call smithEtAl( darkSmith, k(i), SCALE_FACTOR, darkSmithPS )

                write( 6, fmt = '(3e11.3)' ) k(i), lightSmithPS, darkSmithPS

        end do

end program PS
